Heat transfer mechanisms in bubbly Rayleigh-Benard convection 
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The heat transfer mechanism in Rayleigh-Benard convection in a liquid with a mean temperature 
close to its boiling point is studied through numerical simulations with point-like vapor bubbles, 
which are allowed to grow or shrink through evaporation and condensation and which act back on 
the flow both thermally and mechanically. It is shown that the effect of the bubbles is strongly 
dependent on the ratio of the sensible heat to the latent heat as embodied in the Jacob number 
Ja. For very small Ja the bubbles stabilize the flow by absorbing heat in the warmer regions and 
releasing it in the colder regions. With an increase in Ja, the added buoyancy due to the bubble 
growth destabilizes the flow with respect to single-phase convection and considerably increases the 
Nusselt number. 



I. INTRODUCTION 

Thermal convection is an omnipresent phenomenon in 
nature and technology. The idealized version thereof is 
Rayleigh-Benard (RB) convection - a single-phase fluid 
in a closed container heated from below and cooled from 
above. A key question is the dependence of the heat 
transfer rate (as measured by the Nusselt number) for 
given temperature difference between the hot bottom and 
cold top plate (i.e., Rayleigh number), given fluid (i.e., 
Prandtl number), and given aspect ratio. In the last 
two decades there has been tremendous progress on this 
and related questions by experiment, theory and numer- 
ical simulation, see [l| and [2| for a recent review. Most 
of the work focused on RB convection for single-phase 
flow. Various situations in the process and energy in- 
dustries, however, involve convection in the presence of 
phase change, e.g. condensing vapors and boiling liquids. 

The effectiveness of boiling as a heat transfer mecha- 
nism has been known for centuries and the process has 
formed the object of a very large number of studies 3. 
Most of the focus has been on the process by which 
the high thermal resistance opposed by the visco-thermal 
layer adjacent to the hot surface is decreased by the vapor 
bubbles, the two main mechanisms believed to be micro- 
convection and latent heat transport. Another significant 
effect of the bubbles, however, is to promote strong con- 
vective currents in the liquid, thus helping remove the 
heated layer near the hot wall. This aspect of the pro- 
cess forms the object of the present study. 

In an actual experiment all the processes occur at the 
same time and it is next to impossible to separately quan- 
tify their relative importance. Numerical simulation ap- 
pears to be a promising tool for this purpose. Ideally, a 
simulation should be able to resolve individual bubbles 
and follow their evolution but, with the present capabil- 
ities, only so few bubbles can be simulated to this level 
of detail that it would be very difficult to draw conclu- 
sive results [1, 0, @, 0] ■ Therefore one has to fall back on 



point-bubble models in which the interaction of the indi- 
vidual bubbles with the surrounding liquid is parameter- 
ized. This approach has proven valuable in the study of 
turbulence in particle-laden flows (see e.g. d, lioll), in 
liquids with gas - rather than vapor - bubbles [ill . Il2l . [l3| 
and for Taylor-Couette flow with microbubbles inducing 
drag reduction [l^ . 

Many important physical mechanisms have been elu- 
cidated by these means and one may therefore hope that 
similar insights might be achieved by extending this line 
of research accounting for phase change processes, and 
the accompanying bubble growth and collapse, in a simi- 
lar way. Thus, to the fluid-mechanic bubble-liquid inter- 
action model used in our earlier work ([H, EIDj we add 
here models for the heat transfer and phase change be- 
tween the bubbles and the surrounding liquid along the 
lines of Refs. [ll,[3. 

The standard single-phase RB convection under the 
Boussinesq approximation is controlled by the Rayleigh 
number 
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where Th and Tc are the temperatures of the hot (bot- 
tom) and cold (top) plate, respectively, H is the height of 
the convection cylinder, g the gravitational acceleration, 
K the thermal diffusivity of the liquid, v its kinematic vis- 
cosity, and /3 the isobaric thermal expansion coefficient). 
The Prandtl number is defined as 



Pr 



(2) 



and the aspect ratio of the cylinder as the ratio of the 
diameter to the height. In this paper we consider convec- 
tion for which, without bubbles, Ra = 2 x 10"'' and Pr 
= 1.75 (water at 100°C); the aspect ratio is 1/2 and the 
cell is cylindrical. With these parameter values, in the 
absence of bubbles, there is a convection roll with fluid 
rising along one side of the cell and descending along the 
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FIG. 1: Vertical velocity in the plane of symmetry of the 
full cylinder in the absence of bubbles; Ra = 2 x 10^, Pr — 
1.75, Nu = 4.75. As throughout the paper, the velocity is 
made dimensionless by using the free-fall velocity {PgH{Th — 



opposite side (see figure [IJ; the Nusselt number has the 
value 4.75. 

Vapor bubbles introduce a crucial new parameter, the 
Jacob number 



Ja 
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in which L is the latent heat, pv and p the vapor and 
liquid density, respectively, Cp the liquid specific heat and 
Tsat the saturation temperature of the liquid. With the 
parameter values used in this study, hydrostatic pressure 
variations are not sufficient to cause a significant change 
of Tgat, which therefore is taken as a constant equal to the 
average of the hot and cold plate temperatures. Physi- 
cally, J a represents the ratio of the sensible heat to the 
latent heat. A very small Jacob number may be thought 
of as a very large value of the latent heat, which will 
tend to limit the volume change of the bubbles due to 
evaporation or condensation. 

For Ja = the latent heat is effectively infinite and 
bubbles cannot grow or shrink; they maintain their ini- 
tial diameter at nucleation, which we take to be 25 ^m. 
Another control parameter in our model is the total num- 
ber Nij of bubbles in the cylinder. Though in real systems 
this number will fluctuate in time somewhat, here we take 
it as constant: Whenever a bubble reaches the top of the 



cylinder and is removed, a new bubble of the standard 
initial size (25 ^m) is nucleated at the bottom plate at 
some random position. 



II. MODEL 

We study the problem in the standard Boussinesq ap- 
proximation augmented by the momentum and energy 
effects of the bubbles, treated as points. When the vol- 
ume occupied by the bubbles is very small, the liquid 
continuity equation retains the standard incompressible 
form 







(4) 



in which u is the liquid velocity field. The momentum 
equations is 
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where D/Dt is the convcctive derivative, p and T are 
the pressure and temperature, and p — vp the dynamic 
viscosity. The effect of the bubbles has been approxi- 
mated in the standard way by modeling them as point- 
like sources of momentum, which is adequate when the 
volume fraction is small and the bubble radius is smaller 
than the fluid length scales (see e.g. [Tsj). 

The position of the ji-tli bubble is denoted by x„ and 
the force L that it applies on the liquid is modelled as 
(see e.g. 
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in which Rbn is the radius of the n-th bubble and the liq- 
uid acceleration is evaluated at the position of the bubble. 
A similar term multiplied by the vapor, rather than the 
liquid, density has been neglected here. 
The liquid energy equation takes the form 

DT 

pcp— = fcV^T + J2 Qn5{^ - x„) (7) 



where k = npcp is the liquid thermal conductivity and Qn 
is the energy source or sink due to phase change of the n- 
th bubble. We model the thermal exchange between the 
?i-th bubble and the liquid by means of a heat transfer 
coefficient h^n and write 



(8) 



where Tn ~ T{xn,t) is the liquid temperature evaluated 
at the position of the n-th bubble. In writing this relation 
we have used the fact that, for moderate temperature 
differences, phase change is slow and the bubble surface 
remains essentially at saturated conditions (see e.g. [isj). 
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The expressions ([6]) and ([8|) and the use of point sources 
of momentum in ([5|) and of energy in ([7]) assume that 
the bubbles interact only through the average fields but 
not directly, which is a reasonable approximation at the 
vapor volume fractions considered here (see e.g. [13, [3) ■ 

Part of the system energy is carried by the bubble 
phase. If El, denotes the energy of a single bubble, n 
the bubble number density and v the bubble velocity, 
conservation of this component of the system energy is 
expressed by 

^ (nEb) + V • (nEbv) = - ^ Q„5(x - x„) (9) 

n 

where the small pdV contribution has been neglected. 
Adding ([7]) and ([9]) gives an equation for the balance of 
the total system energy, namely 

^ [pCpiT~T,at)+nEb]+V■[pCp{T-T,at)u+nEb^^)] = kV^T. 

(10) 

With the neglect of the vapor mass, the equation of 
motion for each bubble balances added mass, lift, and 
buoyancy, 



in which py = Pv{Tsat)- 

In order to complete the mathematical formulation of 
the problem, definite choices must be made for several 
quantities. Since our bubbles arc small and therefore will 
not deform very much, we take Ca — 1/2, the standard 
potential- flow value for a sphere (see e.g. |24{). indepen- 
dent of the Reynolds number and of non-uniformities of 
the flow [25„ 26, 21,, M- The inviscid calculation of ^ 
gives the same value for the lift coefficient; this value 
appears to be a reasonable estimate even at low to mod- 
erate Reynolds number (see figure 17 of 1231 ) . We model 
the drag coefficient as suggested by [2l|, |29| . 
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in which Ca, Cl, and Cd are the added mass, lift 
and drag coefficients, respectively. The uncertainty with 
which many of the terms of this equation arc known is 
well appreciated in the literature (see e.g. [l^ or our own 
work [20|). Moreover, due to interaction with the wake, 
there rnight be history forces which have been neglected 
in (fTT|) [2l|, [23, . Nevertheless, as written, the equation 
captures the basic effects of drag, buoyancy, and added 
mass which dominate the bubble-liquid interaction. Af- 
ter some rearrangement, the equation becomes 
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The bubble radius Rb is calculated by balancing the la- 
tent heat associated to evaporation or condensation with 
the heat exchanged with the liquid 
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Since the bubble is assumed to be at saturation, pv is a 
constant and this equation can be simplified to the form 

dRb _ hb 
dt Lpv 



{T - T,at) 



(14) 



in which Reb — 2i?f,|v — u|/i^ is the bubble Reynolds 
number. 

We express the heat transfer coefficient hb in terms of 
a single-bubble Nusselt number 
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The dependence of Nub on the parameters of the problem 
is complicated and has been studied by several authors 
(see e.g. MM)- 

In order to make progress we are forced 
to introduce some simplifications. The analysis of [Tsj 
shows that, as a function of the Pcclet number 
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there are essentially two regimes. At low Pcb, Nub is 
approximately independent of Pcb and only depends on 
the Jacob number Wc call this value Nub,o The 

functional relationship Nubfi{ Ja) in this regime has been 
variously parameterized by different authors. Reference 
proposes a general form 
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For the function f{Ja) Ref. [30] (corroborated by the 
more recent results of Ref. [l^ ) proposes 



with which (flSl) becomes 



(671^)1/3 1 3 
16 Ja2/3 + 4 



Nubn 



^ feJaV^^ 12 
2+ +—Ja 

\ IT J IT 



(19) 



(20) 



For very large Peclet numbers, heat transfer is dominated 
by convection and the result is [3l[ 
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FIG. 2: The interpolation (|22|) for the dependence of the 
single-bubble Nusselt number on the Peclet number. 



We combine these two asymptotic forms in a way that 
smoothly interpolates between them: 
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where 



Refs. 31 



~ 2.65 is determined by fitting the results of 
and and the crossover Peclet number Pcc, 
defined by Nub^oo = Nubfi, is Pcc = ttNuI^/A. The 
relation is shown as a function of Peb for Ja ~ 1 
and 10 in figure [H These results can be compared with 
the corresponding ones presented in figure 3 of [l3| and 
are seen to provide an accurate representation of them. 



III. NUSSELT NUMBER 

If the total energy equation ((TO)) is averaged over time 
and integrated over the cylinder volume we find 

{nEbV3 - kd3T)A.t\^^H ^ {nEbV^ - kdsT) A.t\^^Q 

(23) 

where the subscript 3 denotes the vertical direction and 
(. . .)A,t the time and area average. In deriving this re- 
lation we have used the no-slip condition for the liquid 
phase and the assumed adiabaticity of the lateral walls. 
The bubble velocity on the bottom and top plates must 
account for the injection and removal of bubbles and 
therefore cannot be taken to vanish. A similar treatment 
of the bubble energy equation © gives 

{nEbV3)Aj\^^H - {nEbV3)Aj\^^„ = (^Qn^ 

(24) 

where R is the radius of the cylinder. The summation 
in the last term is over all the bubbles contained in the 
system and the average is over time. Using (P^ . can 



equivalently be written as 
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which expresses the obvious fact that any difference be- 
tween the heat conducted out of the bottom plate and 
into the top plate is due to the energy stored in the bub- 
bles. 

In single-phase natural convection the conventional 
definition of the Nusselt numbers Nuc and Nuh at the 
hot and cold plates is 

N^c.h = {d3T)AA.=H...=o (26) 

In the single-phase case this quantity may be considered 
as a total dimensionless heat flux, but this interpretation 
would be incorrect here as it disregards the effect of the 
bubbles. Here the proper quantity to be regarded as the 
total dimensionless heat flux would be 
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as expected. However, since the point of this paper is 
to show the impact of the bubbles on what would be 
considered the heat flux in single-phase convection, it is 
preferable to present our results in terms of Nuh c rather 
than Nl^. 

The definitions (1261) lead to 
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Separate expressions for Nuc and Nuh can be found by 
using another relation which can be derived by multiply- 
ing ([7]) by z — iiJ and integrating over the volume of the 
cylinder with the result 

1 s H 
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in which (. . .)v.t denotes a time and volume average; in 
the following we refer to Nu as the average Nusselt num- 
ber. By using this relation and ((29|) we have 



Nuh = 1- 
and 
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The dimensionlcss heat fluxes ^ can be recon- 
structed by noting that, since bubbles are injected with 
a small velocity and a small radius, the first term in the 
right-hand side of ([M]) is much smaller than the second 
one and therefore, approximately. 



{nEbV3)A,t\ 




Thus, at the hot plate, 

K ~ Nuh (34) 

and, at the cold plate, 

(35) 

in the last step of which use has been made of ([29|) . 

Just as the Nussclt number, the expressions for the 
kinetic and thermal dissipations e„ and eg of standard 
single-phase natural convection are also affected by the 
bubble contribution to the liquid energy equation. These 
modified expressions are derived in the Appendix. 

IV. NUMERICAL METHODS 

Equations ^ and ([7|) have been written in cylindri- 
cal coordinates and discretized using staggered second- 
order-accurate finite difference schemes. The resulting 
algebraic system is solved by a fractional step method 
with the advective terms treated explicitly and the vis- 
cous terms computed implicitly by an approximate fac- 
torization technique (see [l^] for details). The Poisson 
equation that enforces the flow incompressibility is solved 
by a direct procedure which relies on trigonometric ex- 
pansions in the azimuthal direction and the FISHPACK 
package [33| for the radial and axial directions for which, 
therefore, a non-uniform mesh distribution can be used. 
The grid is non-uniform in the radial and axial directions 
and clustered towards the boundaries to adequately re- 
solve the viscous and thermal layers. Following Verzicco 
and Camussi , we used a grid with 33 x 25 x 65 points, 
respectively, in the azimuthal, radial and axial directions 
after having verified that this resolution is sufficient for 
the present Rayleigh and Prandtl numbers. 

Although the code can handle high-order multistep 
schemes, the time advancement of the solution has been 
carried out by a simple second-order Adams-Bashforth 
procedure. For this problem, the most severe limitation 
on the time step size is imposed by the bubble relaxation 
time which, especially for the smallest bubbles, is much 
more stringent than the flow stability condition. 

The only relevant change with respect to the method 
described in [U is the presence of bubble-induced mo- 
mentum and thermal forcings in the governing equations. 
The forcing due to each bubble is located at its center 



and therefore, when ^ and ([7]) are discretized, it has to 
be replaced with an equivalent system of forcings at the 
grid nodes. For this purpose, since in a staggered grid 
arrangement the momentum cells in the three directions 
are all different, the force ([6]) exerted by the bubble is 
first decomposed into its radial, azimuthal and vertical 
components. Each one of these components is then dis- 
tributed by suitable weighing among the 8 vertices of the 
surrounding momentum cell in the same direction. For 
example, for a radial force component / at a position 
n + ^Ar, Oj + 7]Ae, Zk + (Az, with Ar, AO and Az the 
grid spacings and < ^,77, C < 1, the portion attributed 
to the node (r^, Zk) is 

/(1-0(1-^)(1-C)- (36) 

The system of 8 forces thus obtained produces the same 
net resultant and couple as the original bubble force. The 
same strategy has been used for the thermal forcing so 
that the total amount of heat that each bubble exchanges 
with the liquid is preserved. 

The bubble trajectory is computed using the Adams- 
Bashforth scheme for position and the Crank-Nicholson 
scheme for velocity. This latter implicit scheme avoids 
the numerical instability induced by the fast dynamics 
of the smallest bubbles. Equation for the bubble 
radius is integrated explicitly. 

The numerical solver has been validated by monitoring 
the temporal evolution of a single bubble in a quiescent 
flow without a thermal field. Furthermore, our results 
have been compared with the theoretical prediction of the 
lateral force induced on a spherical bubble rising with a 
constant velocity in a viscous fluid near a vertical cylin- 
drical wall. We followed the theoretical method of ref. 
[35! using the free-slip boundary condition for the bub- 
ble surface instead of the no-slip condition used for a rigid 
particle. Another test of the numerical method and its 
implementation is offered by a comparison of the numer- 
ical results for the two sides of Such a comparison 
is shown for a typical case in Fig. [S] 

V. IMPLEMENTATION 

From the numerical point of view, a significant practi- 
cal difficulty of the present problem is the large difference 
between the flow time scale and the times over which 
bubbles grow and collapse. In order to have reasonable 
execution times of our computer code it has been neces- 
sary to limit this difference by adopting a small cylinder 
size; we have taken a height H ~ 17.9 mm and a di- 
ameter 2R ~ 8.94 mm. Furthermore, in order to limit 
the number of spatial cells necessary to resolve the flow 
it is necessary to limit the Rayleigh number, which can 
be achieved by taking a small temperature difference; we 
take Tfi — Tc = 0.25 K. With these values and the physi- 
cal properties of water at 373 K, we have Ra = 2 x 10^. 
Since T^at = ^{Th + TJ, the hot plate is 0.125 K hotter 
than the saturation temperature, which in reality would 
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FIG. 3: The line shows the numerical results for the left- 
hand side of (|29p and the points those for the right-hand side. 
Equality of these two quantities is a stringent of the accuracy 
of the computation. Ni, — 5000. 



not be a superheat sufRcicnt to nucleate bubbles. This is 
another respect in which our model deviates from real- 
ity. On the other hand, since our focus here is the bubble 
effect on the thermal convection, rather than the actual 
heat removal from the plate due to bubble formation, the 
compromise that is forced on us is less damaging than it 
might be in a study of boiling heat transfer. 

The calculation is started without bubbles and run un- 
til the steady state shown in figure [T]is reached. At this 
point 25 /.tm-diameter bubbles are introduced randomly 
throughout the volume of the cylinder attributing to each 
one the local liquid velocity. From this point on, when- 
ever a bubble reaches the top plate, it is removed and 
a new 25^m-diameter bubble is introduced at a random 
position on the bottom plate. The new bubble is placed 
at a height above the plate equal to its radius and it is 
given the local liquid velocity. Bubbles reaching the lat- 
eral vertical wall of the cylinder are assumed to bounce 
elastically. 

In order to avoid possible numerical problems due to 
the disappearance or excessive growth of bubbles, we 
have imposed artificial limits on the minimum and maxi- 
mum bubble diameters equal to 0.82 ^m and 258 /im re- 
spectively. We found however that these limits are never 
approached in our simulations. Since bubbles never con- 
dense completely, the total number of bubbles is constant 
in time. 



VI. RESULTS 

Since bubbles tend to grow in volume in hotter liquid 
region, thus aiding buoyancy, and to condense in colder 
regions, they have a de-stabilizing effect on natural con- 
vection. These effects are clearly the stronger the larger 
the volume change. As explained before, in the present 
model this feature can be controlled by controlling the 
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FIG. 4: Averaged Nusselt number Nu vs Jakob number for 
three different numbers of bubbles. 



Jacob number A very small Jacob number may be 
thought of as a very large latent heat, which will tend to 
limit the volume change of the bubble, while, conversely, 
a large Jacob number would enhance the destabilizing 
effect. 

While this is the major effect, there are other minor 
ones which operate in the opposite direction. For exam- 
ple, bubbles in a hot liquid region, for which T > Tsat, 
will tend to cool the liquid by absorbing heat, and con- 
versely in a colder liquid region. If Ja is very small so 
that the bubble is prevented from growing appreciably, 
this process tends to eliminate the very temperature dif- 
ferences which drive the natural convection in the first 
place. All other things being equal, the break-even point 
between increased buoyancy due to bubble expansion 
and decreased liquid buoyancy due to the bubble-induced 
cooling will be for that value of the Jacob number at 
which the thermal expansion of the bubble equals the 
added weight of the liquid due to the increased density. 
It will be seen from our results that this balance occurs 
for very small J a so that, in most practical situations, 
the balance will tip in favor of the enhanced buoyancy 
effect. 

Figure 2] shows the effect on the average Nusselt num- 
ber Nu = ^{Nuh + Nuc), defined in ([3D]), of adding 
1,000, 5,000 and 10,000 bubbles to the basic single-phase 
RB flow; here, as in all the results shown, the Rayleigh 
number is Ra = 2 x 10^ and Pr ~ 1.75. Figure [S] shows 
the fraction of the bubble contribution 

to the average Nusselt number Nu. The remaining frac- 
tion of the Nusselt number is accounted for by conduc- 
tion and pure convection, i.e. the first two terms in the 
right-hand sides of ([31]) and ([5^ . In both figures the 
horizontal axis is the Jacob number, which we use as a 
control parameter to investigate the effect of the added 
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FIG. 5; Ratio of the bubble source term (|37|l to the average 
Nusseh number (|30p . At smaU Ja, the additional bubble- 
induced buoyant forcing is the dominant effect, while at large 
Ja the bubbles act as direct carriers of heat from the bottom 
to the top. 




FIG. 6: Comparison between the Nusselt number computed 
at the top and at the bottom boundaries for 5,000 bubbles; 
the middle line is the average Nuselt number Nu = | {Nuh + 
Nuc). Here we took Nt = 5000. 



bubble buoyancy. 

For Ja — the bubbles maintain their initial diameter 
at injection at the plate (25 ^m) but, because they are 
kept at Tsat , they cool the hotter liquid regions and heat 
up the cooler ones. As noted before, this behavior tends 
to stabilize the RB convection and is responsible for the 
fact that, while in the absence of bubbles the flow con- 
sists of an annular roll with an approximately horizontal 
axis (Fig.dJ), the addition of Ja = bubbles changes it to 
a toroidal roll with a vertical axis. Because of this stabi- 
lizing effect, the cooling/heating operated by the bubbles 
accounts for a large fraction of the total heat transported 
and, indeed, it can be seen from Fig. \5\ that the bubble 
contribution ([HT]) is very large, up to about 90% of the 
total for the 10,000 bubble case. 



FIG. 7: Average void fraction in the up-flow and in the down- 
flow regions for Nt = 5000 bubbles. 



As Ja is increased, the Nusselt number increases very 
rapidly at first (Fig. |4|) due to the increased convection 
caused by buoyancy. As a consequence, the fraction of 
the total Nusselt number due to the bubbles (Fig. [5]) un- 
dergoes a steep decline. With further increases of Ja, the 
Nusselt number keeps growing but at a more moderate 
rate. The minimum around Ja ~ 0.1 observed in figure 
[5] is due to a change of the flow structure as described 
later. 

Figure [H] shows the Nusselt numbers computed at the 
top and bottom of the cylinder and their average for 5,000 
bubbles; the behavior for the other bubble numbers is 
very similar. As shown by ((29|) . the difference Nuc — Nuh 
is due to the heat exchanged between the bubbles and 
the liquid. As the Jacob number begins to increase, the 
energy absorbed by each bubble per unit time increases 
because of a direct increase in the heat transfer coefficient 
of each individual bubble (see Eq. [201), and an increase in 
the convective component of the bubble heat fiux caused 
by the faster rise velocity of a larger bubble (Eq. [2T|l . 
The moderation in the rate of growth of Nu at larger 
Ja is probably due to the increasing bubble rise velocity 
which limits their residence time in the cylinder. 

By calculating the volume of bubbles located in re- 
gions of positive and negative vertical liquid velocities 
we can look in detail at the effect of the increased buoy- 
ancy. Figure [7] shows the time- and volume-averaged va- 
por volume fractions for 5,000 bubbles as a function of 
the Jacob number. The results for the other cases are 
similar, with smaller void fractions for 1,000 bubbles (for 
J a = 0.35, approximately 0.02% and 0.08%), and larger 
ones for 10,000 bubbles (for Ja ~ 0.35, approximately 
0.16% and 0.36%). It is seen that the void fraction in 
the upflow regions is consistently much larger than in the 
downflow regions, thus providing strong evidence for the 
expected destabilizing effect of the buoyancy provided by 
the bubbles. 

The void fraction reflects the combined effect of bubble 
number and bubble volume and it is interesting to con- 
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FIG. 8: Averaged radius of the bubble computed in the up- 
fiow and down-flow regions for A'^i, = 5000 bubbles. 
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FIG. 9: Averaged bubble numbers in the up-flow and down- 
flow regions for Nt — 5000 bubbles. 



sider these two contributions separately. The volume- 
and time-averaged bubble radius {Rb)v,t, defined by 



1/3 



(i?, 



b/V,t 



bi)t 



(38) 



is shovi^n in Fig. [S] as a function of the Jacob number for 
the case of Fig. [7] vifith 5000 bubbles. As expected, the 
bubble size increases markedly with the Jacob number 
and it tends to be somewhat larger in the hotter liquid 
regions. The time- and volume-averaged fractions of the 
total bubble number in the upflow and downflow regions, 
shown in Fig. [9l indicates a strong tendency for bubbles 
to be in the hotter liquid regions, which is mostly re- 
sponsible for the much larger void fraction in the rising 
liquid. This effect is probably due to fact that the newly 
injected bubbles at the hot plate tend to be swept up into 
the warm liquid by the convection current. 

The results of Fig. [5] for the bubble numbers show that 
the difference between the fractions of bubbles in the up- 
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FIG. 10: Vertical and horizontal cross sections (taken at 
0.05H, 0.5H, and 0.95H, respectively) of the vertical liquid ve- 
locity distribution in the cylinder for Ja — and Nb = 5000 
bubbles. The blue structure near the axis is the descending 
region of the toroidal vortex which prevails for small Jacob 
numbers. The absolute values of the velocities are two or- 
der of magnitude smaller as compared to the two subsequent 
flgures as convection is suppressed at Ja — 0. 



flow and downflow regions is very large for small Jacob 
numbers and tends to decrease as Ja increases. This be- 
havior can be understood looking at the change in the 
flow structure. 

Without bubbles, the cylinder is occupied by a sin- 
gle convective roll which rises along one side and de- 
scends along the opposite side (Fig. [J). A picture of 
the flow for the 5,000 bubbles, Ja = case is shown in 
Fig. [TO] where one vertical and three horizontal cross 
sections color-coded with the vertical velocity field are 
displayed. The blue structure in the proximity of the 
cylinder axis is the descending region of a toroidal vor- 
tex, while the remaining green areas are those where the 
liquid rises, mostly with a smaller velocity, except for a 
few faster zones (yellow and red). It can be seen here 
that the volume occupied by the rising liquid is much 
greater than that occupied by the descending liquid, and 
this circumstance offers a likely explanation of the much 
smaller fraction of bubbles in the latter. 

If the Jacob number is increased to Ja = 0.0935 
(Fig. [TT|) . the toroidal circulation is reinforced with a 
marked increase in the maximum rising and descend- 
ing velocities (note that the color scales in these fig- 
ures are not the same). For a still larger Jacob number, 
Ja = 0.374 (Fig. [T^ the flow has changed back to a 
circulation rising along one side of the cylinder and de- 
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FIG. 13: Fourier modes of the kinetic energy in the angular 
direction for TVj, = 5000 bubbles. Mode corresponds to 
a toroidal vortex and mode 1 to a circulatory motion in a 
vertical region with approximately horizontal axis. 



FIG. 11: Vertical and horizontal cross sections (taken at 
0.05H, 0.5H, and 0.95H, respectively) of the vertical liquid ve- 
locity distribution in the cylinder for Ja = 0.0935 and 5,000 
bubbles. The blue structure near the axis is the descending 
region of the toroidal vortex which prevails for small Jacob 
numbers. 
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FIG. 12: Vertical and horizontal cross sections (taken at 
0.05H, 0.5H, and 0.95H, respectively) of the vertical liquid 
velocity distribution in the cylinder for Ja = 0.371 and 5000 
bubbles. The blue structure near the axis is the descending 
region of the toroidal vortex which prevails for small Jacob 
numbers. 




scending along the opposite one reminiscent of the single- 
phase pattern of Fig. [T] Now the volumes occupied by 
the two streams are more balanced and the difference be- 
tween the number of bubbles in the upflow and downflow 
regions is smaller, as seen in Fig. [91 although the bubble 
fraction in the upflow is still larger than in the downflow. 

These qualitative observations on the flow structure 
can be made quantitative by an analysis of the distribu- 
tion of the liquid kinetic energy among different Fourier 
modes in the angular direction. We define the portion 
En of the kinetic energy pertaining to mode n by 



En = 



rdr 



H 



dz(|u„| 



(39) 



where u„ is the n-th Fourier coefficient (in angular di- 
rection) of the vector velocity field. The mode n = is 
axisymmetric and corresponds to a toroidal circulation 
symmetric about the vertical axis of the cylinder; n = 1 
is a single vortex around an approximately horizontal 
axis, and the higher modes give further information on 
the details of the distribution of the flow over the cross 
section of the cylinder. Results for the ti = 0, 1 and 2 
modes are shown in Figs. [T51 and [HI for 5,000 and 10,000 
bubbles, respectively; the time averaging was carried out 
over the entire duration of the two-phase simulation. The 
values for Ja = are very small, but non-zero. It is seen 
here that, for zero or small Jacob number, most of the 
kinetic energy is in the toroidal mode n = which, for 
5,000 bubbles, reaches a maximum at J a ~ 0.09, which 
is the case shown in Fig. [TlJ For larger values of Ja, 
the energy in the n = mode decreases while that in 
the n = 1 mode rapidly increases giving rise to the flow 
structure exemplified in Fig. [T^l 
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FIG. 14: Fourier modes of the kinetic energy in the angular 
direction for A^i, = 10000 bubbles. Mode corresponds to 
a toroidal vortex and mode 1 to a circulatory motion in a 
vertical region with approximately horizontal axis. 



by the no-slip condition on the cyhnder walls, 

Eu = v{djUidjUi)v = Pg{{T~Tsat)uz)v'\ — |j ^ (f,i • u)^ 

n 

(Al) 

The term ((T — Tsat)u3)v can be eliminated in terms of 
the single-phase Nusselt number at the hot base of the 
cyhnder, given by (pij) . to find 



Ra 
pcpV \ 



n 



(A2) 



in which V = ttR^H is the volume of the cylinder. Al- 
ternatively, in terms of the Nusselt number at the cold 
top of the cylinder. 



VII. SUMMARY AND CONCLUSIONS 

In this paper we have presented a simple model to sim- 
ulate the effect of phase change and two-phase flow on 
natural convection. While, for the reasons given in sec- 
tion |Vl the results must be considered as preliminary, we 
have found that the addition of bubbles has a profound 
effect on the flow structure and on the Nusselt number. 
Bubbles that are prevented from growing by artificially 
maintaining the Jacob number equal to zero (correspond- 
ing to an infinitely large latent heat of vaporization) tend 
to short-circuit temperature non-uniformities and to sta- 
bilize the convective motion. As the Jacob number is 
increased, the added buoyancy due to the bubble growth 
rapidly increases the circulation and the heat transport. 
As the Jacob number is increased further, bubble growth 
is rapid, the residence time short, and the rate of growth 
of the Nusselt number slows down. Correspondingly with 
the increasing Jacob number, the structure of the convec- 
tive flow in the cylinder undergoes significant changes. 



i^^ Ra 
pCpV \ 



(^"c-l) + ^5](f,M-u) 

n 



(A3) 



The thermal dissipation eg is defined in terms of 



e = T --[Th +T,) = T - T,at (A4) 

as eg ~ K(|V0p)y_t. An expression for this quantity may 
be readily obtained by multipling the energy equation by 
9 and averaging over the cylinder volume and time to 
find 



2H 



(A5) 
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APPENDIX A: EXACT RELATIONS FOR THE 
KINETIC AND THERMAL DISSIPATIONS e„ 

AND ee 

Upon multiplying the momentum equation ([5]) by u 
and averaging over the cylinder volume and time, we find. 



where we have used the assumed insulation of the lateral 
walls and the fact that = ^^{Th — Tc) at the bottom 
and top of the cylinder. The temperature gradients can 
be eliminated in terms of the Nusselt numbers Nuh^c to 
find 



kA^ Nuh + Nuc 1 
1P 2 ^ pCpV 



^^(T?i — Tsat)Qn 



which replaces the well-known relation eq 
{kA"^ / H'^)Nu of single-phase RB convection. 
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